/*******************************************************************************
 2HDMC - two-Higgs-doublet model calculator
 Demo program

 http://2hdmc.hepforge.org
*******************************************************************************/
#include "Constraints.h"
#include "DecayTable.h"
#include "HBHS.h"
#include "SM.h"
#include "THDM.h"
#include <iostream>
#include <fstream>

using namespace std;

int main(int argc, char *argv[]) {

	// Reference SM Higgs mass for EW precision observables
	ofstream fout("Br_hgaga_BP1.txt");
	if(!fout){cout<<"can't open file"<<endl;}
	else{fout<<"mh" << "    " << "sba" << "    " << "tb" << "    " << "Gamma(h>bb)" << "    " <<"Br(h>bb)" << "    " <<"Gamma(h>gaga)" << "    " <<"Br(h>gaga)" << "    " <<"Total Gamma(h)" << "    "<< "hb_result" << "    " << "hb_channel"<< "    " << "hb_obsratio" << "    "<< "hs_chisq" << "    " << "hs_pvalue" 	<< "Gamma(H>hh)" << "    " << "Br(H>hh)" << "    " << "Total Gamma(H)" << "    " << "CS_signalH" << "    " << "CS_Hhh" << "    " << "CS_bbgaga" << "    " << endl; }
	
	// Create SM and set parameters
	SM sm;
	sm.set_qmass_pole(6, 172.5);
	sm.set_qmass_pole(5, 4.75);
	sm.set_qmass_pole(4, 1.42);
	sm.set_lmass_pole(3, 1.77684);
	sm.set_alpha(1. / 127.934);
	sm.set_alpha0(1. / 137.0359997);
	sm.set_alpha_s(0.119);
	sm.set_MZ(91.15349);
	sm.set_MW(80.36951);
	sm.set_gamma_Z(2.49581);
	sm.set_gamma_W(2.08856);
	sm.set_GF(1.16637E-5);
	
	// Create 2HDM and set SM parameters
	THDM model;
	model.set_SM(sm);
	
	// Set parameters of the 2HDM in the 'physical' basis
	double mH = 125.;
	double mA = 121.812799;
	double mC = 161.021149;
	//double sba = -0.091551;
	double lambda_6 = 0.;
	double lambda_7 = 0.;
	double m12_2 = 22.299875;
	//double tb = 1.262490;
	
	double gamma_hgaga=0., br_hgaga=0., gamma_total_h=0.;
	double gamma_hbb=0.,   br_hbb=0.  ;
	
	double gamma_Hhh=0, gamma_total_H=0, br_Hhh=0;
	
	double CS_singleH=0, CS_Hhh=0, CS_bbgaga=0; 
	//const HBHSResult *hbhsres_ptr = nullptr;
	HBHS hbhs{};

	for(double loop_tb=4; loop_tb<60; loop_tb=loop_tb+5){
		for(double loop_sba=-0.2; loop_sba<0; loop_sba=loop_sba+0.05){
			for(double loop_mh=20; loop_mh<60; loop_mh=loop_mh+5){
				bool pset = model.set_param_phys(loop_mh, mH, mA, mC, loop_sba, lambda_6, lambda_7, m12_2, loop_tb);
			
				if (!pset) {
					cerr << "The specified parameters are not valid\n";
					return -1;
			    }
			
			    // Set Yukawa couplings to type II
			    model.set_yukawas_type(3);
			
			   // Prepare to calculate observables
			    Constraints constr(model);
			
			    double S, T, U, V, W, X;
			
			    constr.oblique_param(mH, S, T, U, V, W, X);
			
			    if(!constr.check_stability()){ 
				//	cout << " Not pass potential stability:  mh "<<loop_mh<<" tb " << loop_tb <<" sinbma " <<loop_sba << endl; 
					continue; 
				}
			    if(!constr.check_unitarity()){ 
				//	cout << " Not pass Tree-level unitarity: mh "<<loop_mh<<" tb " << loop_tb <<" sinbma " <<loop_sba << endl; 
					continue; 
				}
			    if(!constr.check_perturbativity()){ 
				//	cout << " Not pass Perturbativity:  mh "<<loop_mh<<" tb " << loop_tb <<" sinbma " <<loop_sba << endl; 
					continue; 
				}
			
			
				const HBHSResult hbhs_result = hbhs.check(model);
				if(hbhs_result.hb.result[0]!=1 || hbhs_result.hb.result[1]!=1 || hbhs_result.hb.result[2]!=1 || hbhs_result.hb.result[3]!=1 ){
					cout << " Not pass HB: mh "<<loop_mh<<" tb " << loop_tb << " sinbma "<< loop_sba << " hb_result " << hbhs_result.hb.result[0] << " " << hbhs_result.hb.result[1] << " " << hbhs_result.hb.result[2] << " " <<  hbhs_result.hb.result[3] << endl;
				    continue;
				}
//				if(hbhs_result.hs.chisq<0.05){
//				//    cout << " Not pass HS: mh "<<loop_mh<<" tb " << loop_tb << " sinbma "<<loop_sba << " pvalue " << hbhs_result.hs.pvalue << endl; 
//				    continue;
//				}
			
				// Show Cross section for 13 TeV
				cout << " pass all: mh "<<loop_mh<<" tb " << loop_tb << " sinbma"<<loop_sba << " pvalue" << hbhs_result.hs.pvalue << endl; 


				// Prepare to calculate decay widths
				DecayTable table(model);
				gamma_hgaga  =table.get_gamma_hgaga(1);
				gamma_hbb    =table.get_gamma_hdd(1,3,3);
				gamma_total_h=table.get_gammatot_h(1);
				br_hgaga     =gamma_hgaga/gamma_total_h;
				br_hbb       =gamma_hbb  /gamma_total_h;

				gamma_Hhh    =table.get_gamma_hhh(2,1,1);
				gamma_total_H=table.get_gammatot_h(2);
				br_Hhh       =gamma_Hhh  /gamma_total_H;
			
                CS_singleH = hbhs_result.hb.CS_LHC13[1].find("singleH")->second;
				CS_Hhh     = CS_singleH * br_Hhh;
				CS_bbgaga  = CS_Hhh * br_hbb * br_hgaga;

				fout << loop_mh << "    " << loop_sba << "    " << loop_tb << "    " <<gamma_hbb << "    " << br_hbb << "    "  << gamma_hgaga<< "    " << br_hgaga<< "    " << gamma_total_h << "      " << hbhs_result.hb.result[0] <<"      " << hbhs_result.hb.channel[0] << "      "  << hbhs_result.hb.obsratio[0] << "      "<< hbhs_result.hs.chisq << "    " << hbhs_result.hs.pvalue << gamma_Hhh << "    " << br_Hhh << "    " << gamma_total_H << "    " << CS_singleH << "    " << CS_Hhh << "    " << CS_bbgaga << "    " << endl; 
				// Write output to LesHouches file
				//model.write_LesHouches("test_out.lha", 1, 0, 1, hbhsres_ptr);
			}
		}
	}
}
